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that no spurious states appear in the discretization process. 
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I. INTRODUCTION 



There has been a strong interest in the solution of the time-independent Dirac equation 
in the last few decades motivated mostly by new advances in molecular and nuclear physics. 
More recently, the fields of laser-matter interaction and optics have also considered this 
equation because new developments have led to experimental facilities reaching laser inten- 
sities above 10 20 W-cm -2 jl]. The mathematical description of an electron subjected to such 
intense electromagnetic fields necessitates a relativistic treatment {2 — 4 1 and thus, theoretical 
efforts should be based on the Dirac equation instead of the non-relativistic Schrodinger 
equation. 

It is well-known that finding a solution of the Dirac equation is very challenging because 
it has an intricate matrix structure. This complicates analytical approaches and closed-form 
solutions can be found only for highly symmetric systems. For this reason, a numerical 
treatment is required to study more realistic physical processes occurring in molecules or 
heavy ion collisions. However, the development of an accurate numerical approach is also 
difficult because the Dirac spectrum is not bounded from below (and above). This precludes 
the "naive" generalization of well-known methods used to solve the Schrodinger equation 
such as the Galerkin or Rayleigh-Ritz methods. These are based on minimization principles 
and thus, work rigorously only if the spectrum of the differential operator has a lower (or 
upper) bound. Trying to solve the Dirac equation with these methods can thus be very 
dangerous as the Dirac spectrum can be altered by negative energy contributions. This 
problem is named the "variational collapse" and leads to the appearance of spurious states 
in the approximated spectrum. The spurious states are eigenvalues which do not belong 
to the spectrum of the continuous operator and which appear in the discretization process. 
More precisely, let A be an eigenvalue of the Dirac operator H in the mass gap (—mc 2 , mc 2 ) 
(corresponding to bound states) and be its point spectrum in the mass gap, that is the set 
of all A's. Numerically, we approximate the operator H by H n such that lim n ^ 00 H n = H 
(here, n is the dimension of the subspace on which the operator is projected or, loosely 
speaking, the dimension of the Dirac operator matrix once the problem is discretized). The 
discretized Dirac operator H n has eigenvalues given by A n 6 a ^ . The set of spurious states 
a s ~ C (Jfr is defined by the set of all eigenvalues in av> for which lim^oo A n ^ a p. 

There have been many successful attempts in the literature to circumvent variational 



lapse by adapting minimization techniques and basis set expansions to the Dirac operator 



-Il4j. Usually, these techniques can be classified into one of three main categories 12]: 



modification of basis functions, utilization of an operator that has a lower bound but the 
same spectrum as the Dirac operator and transformation of the Dirac operator. The methods 
utilized in this work fall into the first and third categories; they are the Rayleigh-Ritz method 
with kinematically balanced basis functions fill . \l£ 



on a min-max principle 



17 



18| 



16] and the variational method based 



The discretization of the Dirac equation in these two cases proceeds by the utilization 
of a basis set expansion which allows, in principle, a very good accuracy. These techniques 
have been exploited extensively in the non-relativistic case to solve the time-independent 
Schrodinger equation 19] and allow making accurate predictions for the non-relativistic 
spectra of molecules and for other observables. In this work, a B-spline set of basis functions 
will be utilized. This choice is motivated by some interesting properties of B-splines: they 
have compact support (leading to sparse matrix structures), they are very flexible in terms 
of element size and continuity conditions (both are determined by their knot vector), they 
are positive definite and finally, they are linearly independent (they form a complete basis). 



For these reasons, B-splines have also been widely applied to the non-relativistic (see 20 1 

24 1 



18 



21 



for a review on the use of B-spline in molecular physics) and to the relativistic 
cases. Here, we combine B-spline basis functions with the two numerical schemes described 
previously. 

We use these methods to investigate two particular systems: the diatomic molecule H^ 
and the quasi-molecule Tli2 79+ . This is accomplished by computing the relativistic spec- 
trum of the two-centre Coulomb problem in the fixed nuclei approximation. The rationale 
for studying these systems is twofold. First, they are physically relevant in many fields of 
physics. The molecule H^ is very important in chemical physics and its relativistic correc- 



tions, albeit very small, have been the subject of many studies 25M29] . On the other hand, 



the quasi-molecule Tli2 79+ is not stable and dissociates rapidly. However, it is pertinent 
in heavy ion collisions at intermediate energy, where processes such as charge transfer and 
electron-positron pair production are investigated 30|, [3l] , and in high intensity laser-matter 
interaction, where pair production and Quantum Electrodynamics (QED) processes could 



be enhanced in the presence of heavy nuclei 



32|. 



The second reason to look at these systems is that their ground state eigenvalue has 



already been computed and thus, the numerical results obtained from our analysis can 



be compared to results 



different analytica 



24 



26 



28 



29 



35 



rom the literature. More precisely, they have been studied using 



34J and numerical approaches such as the Rayleigh-Ritz scheme 



36], the variational scheme based on the min-max principle [18] and 



finite difference methods 



25l . 127] . Very accurate results for the ground state of diatomic 



molecules were obtained in these analyses. However, the whole spectrum is rarely discussed 
(an exception to this is found in 0, [|8]) and some of these methods (especially those based 
on the "naive" Rayleigh-Ritz method) could potentially lead to the appearance of spurious 
states. For instance, this problem was discussed in 24J and a technique for identifying these 
artifacts was described. However, this can be cumbersome when one is interested in sums 
over intermediate states such as those required in radiative QED corrections. For these 
calculations, it is certainly more efficient to have a numerical scheme free from spurious 
states from the outset. 

In this work, we are presenting and comparing two numerical methods that use a B- 
spline basis set expansion to compute the relativistic spectrum of the two-centre problem. 
The first one is the min-max variational method. The second one is the Rayleigh-Ritz 
method combined with kinematically balanced basis function. In Section [Til t ne variational 
formulation of both numerical methods is presented and the choice of basis functions is 
described. The numerical results are displayed in Section III II where some values for the 
spectra of diatomic molecules are shown along with a discussion of spurious states. The 
conclusion is found in Section [TV] 



II. NUMERICAL METHODS 



The Dirac equation describes the relativistic dynamics of spin-| particles (fermions) like 
the electron. In this work, we consider specifically the single particle static Dirac equation 
without vector potential given by 39] 



Hip(x) = Eip(x) with H = col ■ p + mc 2 f3 + V(x)l^ 



(1) 



where H is the Hamiltonian operator, p = — iV is the momentum operator, c is the light 
velocity, m is the electron mass, E is the electron energy and ip G L 2 (R 3 , C 4 ) is a four-spinor. 
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The matrix structure is given by a and (3 which are four by four matrices given by 





<7j 




h o 


OCi = 




and /3 = 








-h 



(2) 



where <7j are the usual Pauli matrices. The latter are 





1 




-I 




1 


cr x = 




, CTy = 




and a z = 






1 


i 




-1 



(3) 



Written in this way, the Dirac equation is in the Dirac representation. This equation gives 
a consistent description of single bound electrons within the fixed nuclei approximation, i.e. 
when the effect of the nuclei is included in the potential term V and the nuclei are fixed 
in space. This is valid when the masses of the nuclei are much larger than the mass of the 
electron, which will always be the case for the systems considered in this study. 

The main goal of this work is to calculate approximate solutions of ([T]). To achieve this, 
it is convenient to write the four-spinor as if)(x) = [4>(x), x{ x )\ where (f>(x) and x{ x ) are 
two bi-spinors called the large and small components, respectively. The Dirac equation then 
becomes 



V(x) + mc 2 
R 



R 




(j){x) 




4>{x) 








= E 




V(x) — mc 2 











(4) 



where we defined R = —iccr ■ V. This last equation is the common starting point for the 
numerical methods that follow. As will be seen later, it is also handy to decompose the 
latter into two coupled equations as 



Rx(x) — [E — mc 2 — V(x)](j)(x), 
R(j)(x) = [E + mc 2 - V{x)]x{x). 



(5) 
(6) 



The small component can then be written in terms of the large component, yielding 

R 



By substitution, we get 

R 



X(x) 



Rct){x) 



E + mc 2 — V(x) 



(j){x) 



(7) 



E + mc 2 — V(x) 



[E -mc 2 - V{x)]<t>(x) 



(8) 
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for the large component, which belongs to G L 2 (IR 3 ,C 2 ). Note here that the latter pro- 
cedure can also be implemented as a Foldy-Wouthuysen transformation 40] . These two 
relations will be important for the analysis that follows. 

For a diatomic molecule, the case considered in this study, the static potential is 

Zi Z2 



V(x) 



(9) 



|x + Rz\ |x — Rz\ 

where Z\^ are the atomic electric charges, 2R is the inter-atomic distance and z is a unit 
vector in the z-coordinate direction. It represents the static Coulomb interaction of two 
point-like nuclei with an electron in the fixed nuclei approximation. This potential is axially 
symmetric so the number of dimensions can be reduced by one: the azimuthal coordinate 
dependence can be treated analytically by factorization. Thus, the four-spinor in cylindrical 



symmetry reads 29 



3 



(10) 



MZ,v)e iUz+1/2)e 
iXi(^v)e i(jz ~ 1/2)d 
^X2(^v)e iUz+1/2)e 

where j z is the angular momentum projection on the z-axis (it can take one of the values 
j z = — |, — |, — |, |, |, |, ...) and 77, £ are prolate spheroidal coordinates. This choice of 
coordinate system is very convenient for the numerical implementation because the Coulomb 
singularities are situated on the domain boundaries. Also, it has already been utilized in 
accurate evaluation of the diatomic ground state energy in the relativistic 28|, |29[ and non- 



relativistic 



42 



44j cases. For these reasons, these coordinates will be used throughout this 



work even though it was argued in |24j that Cassini coordinates could provide slightly more 
accurate results. The prolate spheroidal coordinates are defined as 

(11) 
(12) 

(13) 



x = R [(£ 2 - 1)(1 -rf)Y cos6, 
y = R[(e-l)(l-ri 2 )]hm6, 
z = R£rj, 



where £ G [1, 00), T] G [—1, 1] and 9 = [0, 2n] (azimuthal angle). The Coulomb potential in 
these coordinates becomes 

Z\ Z 2 



(14) 



We can now start discussing the numerical methods utilized in this work to calculate the 
spectrum of diatomic molecules. 



A. Min-max method 



The first method described in this work was developed in 13, [18|, I45M47I] and is a weak 
formulation for operators with gaps in their spectrum. The main idea is to find the critical 
points of the Rayleigh-Ritz coefficient by using a min-max principle. More precisely, it was 



rigorously proven that the sequence of values defined by 



45] 



inf 

G subspace of F + 
dim G = k 



\H\ 



sup ... 

Ve(GeF_)\{o} {tpm 



(15) 



give the actual eigenvalues of H in the mass gap, if certain assumptions are fulfilled (for in- 
stance, the first eigenvalue should obey Ai > -mc 2 and the potential V should be Coulomb- 
like). Here, F + _ are two well-defined orthogonal subspaces of F C L 2 (M 3 ,C 4 ). This result 
is very general and can be applied to any space decompositions F + © F-. For practical 
purpose, it is convenient to use one that splits the large and small components of the Dirac 
equation as [45 ] 





® L 2 (M 3 ,C 2 ^ 



(16) 



In this setting, the maximization is performed exactly by the relation ([71) which relates the 
large and small components 45(. Then, only the minimization remains and this step is 
carried out numerically. 

The latter can be performed by minimizing the energy over any couple (E, <fi) obeying 
the functional equation given by 



A[E,4>] = / d 3 x 



E + mc 2 -V 



[E-mc 2 - V]\</>\' 



0. 



(17) 



This functional equation is obtained from (J7J) and (jSJ) by multiplying the latter by $ on 
the left, by integrating on space (using integration by parts) and by using the divergence 
theorem. Also, the wave function should vanish faster than ~ 3 at infinity. This is the 



case when V is a Coulomb-like potential because the corresponding wave function vanishes 



as 



~ e r when r — > oo 



48|. 



Therefore, the last equation gives a realization of the min-max principle. Moreover, it 
is shown from the preceding procedure that it predicts the same spectrum in the mass 
gap (—mc 2 ,mc 2 ) as that of the Dirac equation 17, is[ 45-47]. This formulation also allows 
discretization schemes that use a basis function discretization without "variational collapse" : 
the calculated energy spectrum is bounded in the mass gap and does not fall into the negative 
energy continuum. Finally, spurious states does not appear in the calculated spectrum 
without adding any additional conditions [17J, making for a very robust numerical method. 
In the following, we describe the discretization of (IT7|) by using a set of basis functions. 



1. Basis set expansion 



The discretization of ( I17|) with the potential in proceeds by expanding the wave 
function over a set of basis functions. Thus, the bi-spinor can be written as 



N 



'1,2 



(18) 



n=l 



(1 2) (1 2) 

where a n ' are the basis expansion coefficients and Bn ' (£,ry) are the basis functions (to 
be defined later), for components 1 and 2, respectively. 

The explicit expression of (|17p in discretized form depends on the potential considered, 
on the coordinate choice and is a complicated functional of basis functions (some examples 
can be found in 

HQ)- Th e equation A[E, <p] = 0, once discretized, generally has the form 



2JV 



where A(E) is now a matrix, Ak(E) is its fc'th eigenvalue and 



(19) 



,(1,2) 
l k,i 



hi 



for i < N 



.(2) 



aii_ N for i > N 



(20) 



its eigenvector. For cylindrically symmetric systems expressed in prolate spheroidal coor- 
dinates, the case considered in this study, A(E) becomes a 2N x 2N matrix having the 
following structure: 

" A n (E) A l2 (E) 
Al(E) A 22 (E) 



ME) 



(21) 



8 



where An, A22 and A12 are N x N matrices for which explicit expressions can be found in 
Appendix lAl 

Therefore, solving the non-linear (A(E) depends on the energy E) eigenvalue problem 
fjTTJ]) gives an approximation of the energy E and eigenfunctions: the energy of the fc'th 
bound state is a solution of Aj.(E) = where A& is the fc'th eigenvalue of A(E) while wave 
function coefficients in the basis expansion are the A(E) matrix eigenvector coefficients. It 
can easily be demonstrated that Ak(E) are monotonically decreasing functions 18] which 
implies that they have only one root, i.e. only one value of E = E TOOt for which Ak(E mot ) = 0. 
This problem can thus be solved by iteration or any other root finding algorithm. 

B. Rayleigh-Ritz method 

The Rayleigh-Ritz method is well-known and has been studied extensively for both the 
relativistic and non-relativistic cases (see HQ for instance). Starting from (jlj), we can 
multiply by ip> on the left and integrate on space to get another functional equation given 
by 

J d 3 x { [mc 2 + V] \4>\ 2 + (cf)\Rx) + {x\R<f>) + [V- mc 2 ]\ X \ 2 } = 

E j d 3 x{|0| 2 +| X | 2 }, (22) 

which is just an explicit way of writing the well-known Rayleigh-Ritz functional equation 
H = (ip\H\tp) / . The notation (-|-) stands for the Hermitian inner product. In the 
following, we define two operators C and S by 

C[^} = Jd 3 x{ [mc 2 + V] \<P\ 2 + (<f>\R X ) + (x\m + [V~ mc 2 ]\ X \ 2 } , (23) 

S[1>] = I rf 3 x{|0| 2 +| X | 2 }. (24) 

A numerical scheme can be developed from these equations by discretizing the wave function 
over a set of basis functions. 

For a bounded operator (like the Schrodinger operator), the best estimate for the eigen- 
pairs is obtained by a minimization procedure because the Rayleigh-Ritz quotient forms an 
upper bound for the eigenenergy E (if the spectrum is bounded from below). Moreover, 
it can be shown that the minimum of the quotient converges towards the exact eigenpair 
as the number of basis function N — > 00 (see and references therein). If the operator 



is not bounded, as in the case of the Dirac equation, the quantity H does not necessarily 
form an upper bound, although it is still a stationary point (as seen in the previous section, 
the eigenvalues can actually be characterized by a min-max principle). For this reason, 
the convergence of this approach is not guaranteed because the stationary point is a saddle 
point, and spurious states may appear. Therefore, a modification of the method is required 
to improve the convergence. The strategy we are using in this paper is the Kinematically 
Balanced Basis functions (KBBF) described in the next section. 

1. Basis set expansion 

The discretization of fl22|) is very similar to the one in the last section and proceeds 
by expanding the wave function over a set of basis functions. In this case, one writes the 
bi-spinors as 

N 

<Pl,2^V) = J2 a n' 2)B n 1 ' 2) (^V), (25) 
n=l 

N 

Xl,2^V) = Y, C n' 2)X n' 2) ^V), (26) 
n=l 

where oin \ Cn ^ are the basis expansion coefficients and Bn' 2 \C,,f])^n' 2 \C,v) are the 
basis functions (to be defined later), for components 1 and 2, respectively. In the naive 
Rayleigh-Ritz method, the basis functions for both spinors are the same, that is X n ' (£, r/) = 
Bn ' (£,,v)- Substituting fl25|) and f )26|) in f )22|) . we obtain a generalized eigenvalue problem 
in the form of 

Ca = £Sa, (27) 

where the generalized eigenvector a = [a£ , an \ af\ an , , ch , c[ \ Cn ] con- 
tains the basis function expansion coefficients. Finding the solution of this last equation 
corresponds to an extremization on the trial function parameters {af~' 2 * > and cf' 2 ^) and yields 
an approximation of the eigenenergies and eigenf unctions. The functionals C[i/j] and S[tp] 

10 



become AN x AN matrices having the following general structure 



^11 





p(3) 
^11 


p(3) I 





r (l) 


p(3) 
°21 


r (3) 
*^22 


p(3)T 


p(3)T 
^21 


^11 





r (3)T 
*^12 


r-(3)T 

*^22 





p(2) 
*^22 . 



(28) 



q(l) 







'22 








o si? 

s 









(2) 
22 



(29) 



where the matrices C^' 2 ' 3 ^ and S^' 2 '* are N x N matrices. Their explicit expressions are 
given in Appendix [B] 

In the KBBF technique, a special choice of basis functions is made: the basis for the large 
and small components are related by a transformation guaranteeing that the non-relativistic 
equation is recovered in the limit c — > oo. This is achieved by considering the relation 
between the small and large component given in (CO). Substituting ( 1251) in the latter, we get 



R 



(30) 



Note here that substituting this equation in the Rayleigh-Ritz coefficient allows recovering 
the min-max method and its discretization scheme described in the last section. Then, the 
eigensolutions have to be calculated by using an iteration procedure because the functional 
has an intrinsic dependence on the eigenenergy. This is circumvented by neglecting the space 
dependence of the potential over the support of each basis function and redefining the basis 
expansion coefficients as 

1 



.(1,2) 



E + mc 2 — V n 
1 



,(1,2) 



,(1,2) 



(31) 
(32) 



E + mc 2 - V(Z,ri) n ' 
where 14 is a constant coefficient representing the contribution of the potential on the 
support of the basis function n. Thus, we obtain 



R 



11 



(33) 



where the factor 1 / 2mc 2 was included for numerical convenience (it also allows to recover the 
non-relativistic limit exactly if the energy is shifted by mc 2 ). In some sense, the eigenenergy 
and space dependence of the prefactor 1/(E + mc 2 — V) is encoded in the coefficients c« . 
No iteration procedure is required as in the min-max method but the eigenvalue problem 
is larger. Also, the relation between the small and large components is only approximate 
because the constants c n ' have neither energy nor space dependence. However, we expect 
this relation to converge toward the exact one as the number of basis functions is increased 
and their support decreases. In that case, neglecting the spatial variation of the potential 
becomes a better approximation. In the limit iV — > oo, we have 

X(t,v) = ^m,ril (34) 

where f is a bi-spinor. This implies that in the KBBF, the extremization of the 
Rayleigh quotient is performed on /(£, if) rather than on x(£, rj) as in the "naive" Rayleigh- 
Ritz method. The former is consistent with the min-max principle exposed previously and 
the stationary point (or Euler-Lagrange equation) is given by 

(»> 

This equation can also be obtained from a unitary transformation of the Dirac equation 



121 ] . Therefore, in the continuous limit, the exact solution is recovered from the min-max 
principle, which establishes that the two approaches are consistent with each other in that 
limit. 

Explicitly, the basis function expansion is given in prolate spheroidal coordinates by (here, 
we dropped the basis function argument for simplicity) 

N 

<Pi,2^v) = J2at 2) B^ 2 \ (36) 

n=l 

N ( n ^1 

XM> 1) = 2^ E ^ [- d r ~ 7] - &d,BS> , (37) 



71=1 

N 



L(2) 
I c n 


-d r - 








r . 


La) 

I L n 


-d r + 








r . 



n=l ^ ) 

where d r and d z are given in ()A4j) and (1A5|) respectively, while r = R [(^ 2 — 1)(1 — rf)] 1 -. 
These last expression were obtained by expressing the operator R explicitly in prolate 
spheroidal coordinates. 
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Substituting (136|) to (1381) in (l22l) . we also obtain a generalized eigenvalue problem in the 
form of (1271). In this case, the matrix structure is 



C 



p(l) 





p(3) 
^11 


r (3) 1 
^12 





r (l) 

^22 


r (3) 
^21 


r (3) 
^22 


n (3)T 
^11 


r (3)T 
^21 


^11 


r (2) 
^12 


r (3)T 
^12 


p(3)T 
*^22 


p(2)T 
^12 


p(2) 
*^22 



(39) 










q(i) 

°22 
















c(2) o(2) 

°11 D 12 

c(2)T «(2) 

°12 °22 



(40) 



and explicit expressions can be found in Appendix The matrix C and S are 4N x 4iV 



matrices while the other components (C^ 1 ' 2,3 '' and S\j' z> ) are N x N matrices. 

The analytical analysis of convergence of this method and the proof that it does not have 
spurious solutions is clearly a non-trivial matter but was discussed in ll) for the L-spinors 
basis functions. In our case, these properties will be verified empirically by looking at the 
numerical results while a careful analysis of the method is currently under investigation. 



:(i,2h 



Basis functions 



Throughout this work, B-spline basis functions are used (a description of these func- 

n 

tions can be found in [20j). This choice is favored over other techniques because it can be 
easily implemented and because B-splines have compact support, leading to sparse matrix 
structures. This allows using powerful numerical routines for the calculation of eigenvalues. 
More important is the fact that B-splines are linearly independent and form a complete 
basis, which is a necessary condition for the convergence of the Rayleigh-Ritz method for 
both eigenenergies and eigenfunctions ( ll] and references therein). It is also an important 
requirement to avoid errors of order 1/c 4 in the Rayleigh-Ritz bounds which may induce 
spurious states in certain circumstances j^]. 

B-splines basis functions have been studied extensively for solving the time-independent 
Dirac equation because of these important properties. However, most of these studies consid- 



ered atoms or atomic-like systems 



18 



21 



23| , although recently, they were used for diatomic 



molecules 



18 



24]. 
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B-splines are fully determined by their order k^ iV and knot vector using the iterative 



formula 



20 



49] 



(*) = * k . b\-\x) + ti+k X fc^s), (41) 

and initial conditions 

b](x) = 1 for U < x < t i+ i and b] = otherwise, (42) 

where ij's are knot coordinates. The number of knots at a given coordinate determines the 
continuity condition at that point. Therefore, the number of knots should be maximal at 
singular points (at the Coulomb singularity position for instance) to allow for a discontinuous 
behavior of the wave function. Throughout this work, the knot vectors are given by the 
sequences 

1 — £l = ••• — 6c 5 < < ••• < £n 5 +l — ••• = £n e +fc c — (max, (43) 

-1 = rii = ... = r) kn < T] kv+1 < ... < r) nv+1 = ... = r] nv+kv = 1. (44) 

Here, n^ )V are the number of spline functions in £ and 77 coordinates, respectively. The 
knot coordinates can be chosen arbitrarily in the domain under consideration. However, to 
improve accuracy, an exponential sequence with smaller intervals close to the singularities 
is used in this study. The knot sequences and domain structure for diatomic molecules are 
depicted in Figure [H 

The basis function can then be written as the tensor product of B-spline functions as 

B^ 2) ^v) = G^ 2 \^ V )b^(0b-\v), (45) 



where n = G Z 2 , % e [l,n^] and j G [1,^]. The overall factor is defined as [28l. 129] 

G (i ' 2) ^v) = [(e-m-v 2 )]^, (46) 

where 

/ii = 3z~ \ and /i 2 = j z + ^. (47) 

This factor accounts for the angular momentum dependence (remember that j z is the angular 
momentum projection on the z-axis). Moreover, it allows having well-defined integrals in 
the functionals, allowing a better convergence of the method. 
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n 

a 



1 ■- HQ: 



-1 | H, 6^ 



FIG. 1. Discretization of the domain. A circle represents a point where there is a Coulomb 
singularity. The domain under consideration is discretized into a certain number of elements while 
each element is subdivided into smaller regions for the space integration. Note that we are using 
an exponential size distribution with smaller elements close to the Coulomb singularities. There is 
a knot point at every intersection of two dotted lines. 



D. Details of the calculation 



The construction of the matrices appearing in the two numerical methods involves several 
integrals extending over the whole domain. However, because B-splines are compact, the 
integration domains are reduced to the support of each basis function, which are regions 
having x k v elements or less. The integrals are evaluated numerically using the Gauss- 
Legendre quadrature rule. 

The boundary conditions are chosen as 0(£ m ax, v) = an d x(£max,^) = 0. Using B- 
splines, this can be implemented easily by setting & n? (0 = and by considering only — 1 
B-spline functions in £ coordinates. The other boundaries are free. Strictly speaking, these 
boundary conditions may lead to some numerical problems because of the occurrence of the 



Klein paradox 21]: they correspond to the confining of the electron in a box surrounded 



by an infinite potential barrier at £ = £ max . However, similar conditions have been used 
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successfully in 



28 



29 



36[ to obtain very accurate results. Also, it was argued in [18|, [22 1 



that using these zero boundary conditions have negligible effects on the solution if the domain 
is large enough. 

The boundary conditions on the nuclei are more subtle. It is argued in 50] that B-spline 



basis expansion are plagued with non-physical solutions related to a wrong treatment of 
these boundary conditions. A remedy to this problem is also proposed but it is not clear to 
us how this technique can be applied to the diatomic molecule case. However, it is verified a 
posteriori by looking at the spectrum that our boundary conditions does not induce spurious 
states. 

The code performing the calculation is parallelized by using the domain decomposition 
strategy described in 51] . For better performance, the ScaLAPACK library is utilized to 
solve the eigenvalue and generalized eigenvalue problems. In the Rayleigh-Ritz method, 
the latter yields the whole energy spectrum and eigenfunction in one calculation. For the 
min-max method, only one eigenenergy can be calculated at a time because each evaluation 
necessitates a solution of Ak(E) = 0. The latter is solved by a root-finding algorithm based 



on Brent's method 



52]. 



III. RESULTS AND DISCUSSION 

In this section, the results obtained for both numerical methods are presented. First, 
the convergence of the method is analyzed. Then, the spectra of diatomic molecules are 
presented and finally, the absence of spurious states is discussed. 



A. Convergence of the method 

In this section, we are investigating the convergence of our numerical methods. More 
specifically, we study and calculate the ground state of dithorium (Th^ 79-1- which has Z\p = 
90) and dihydrogen (H^" which has Z 12 = 1)). The semi inter atomic distance is set to 
R = ^ w 0.011111 a.u. for dithorium and to R ~ 1.000 a.u. for dihydrogen while the 
angular momentum is taken as j z = 1/2. The results for the calculation of the ground state 
binding energy using B-splines of order 7 and different mesh sizes are shown in Table [J and 
[Til for and Th^ 79 " 1 ", respectively. 
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The results presented in this table show the convergence of the method as the number 
of elements is increased. The results obtained are very accurate, although there is a small 
relative difference d r between our results and the results presented in j^. For Hjj~, we 
obtain d r ~ 2.2 x 10 _8 % and d r ~ 1.9 x 10~ 8 % for the min-max and Rayleigh-Ritz methods, 
respectively, while for Tli2 79+ , we have d r ps 6.1 x 10~ 4 % and d r fa 0.97 x 10~ 4 % for the 
min-max and Rayleigh-Ritz methods, respectively. These differences can be explained by a 
different choice of boundary conditions, different element formulation and different treatment 
of the Coulomb singularity. 

In all cases, the convergence is from above suggesting that both methods are consistent 
with the min-max principle. The convergence in the case of dithorium however is much 
slower than for dihydrogen. One possible reason explaining this discrepancy is the behavior 
of B-splines close to the Coulomb singularities. In that region, the wave function should 
behave like if) ~ r 1 ^ +l1 ' 2 (obtained in atomic calculations) where 




71,2 = A / + - -«% 2 2 . (48) 



and r 12 are the distances from nuclei 1 and 2. In ground state calculations, we have j z = 1/2 
and thus, < 7^2 < 1 for Z% t 2 < 137. Therefore, the wave function has a non-integer power- 
law behavior close to the singularity. The B-spline basis functions, being polynomial with 
integer powers, are unable to reproduce exactly this feature. Moreover, we have that 

7h « 0.999947 and ip ~ r -°-°°««» (49) 
7Th ~ 0.568664 and ip ~ r^' 431336 , ( 5Q ) 

where 7H,Th are the gamma associated with a hydrogen and thorium atom. It is clear from 
this that the behavior of the wave function is much closer to a power law for dihydrogen 
and therefore, is better reproduced by the B-splines and also, has a faster convergence. 

One possible cure to this is to use another prefactor in the basis function that captures 
the correct behavior. For instance, it was proposed to multiply the basis functions in ( 1451) 



by |28|,|29|,|36| 



G ' (e ^ )=r -l+Ti r -l+7 2; (51) 

with 

r 1 = ^ + V )R, r 2 = {t-rj)R. (52) 
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The main issue with this technique is that derivatives in the functionals become singular at 
the nuclei position. To cope with this, a singular coordinate transformation can be performed 
that allows transforming the singular non-integer behavior near the nuclei to a polynomial 



approximation 



28 



361 ] . The latter can then be fitted more accurately with a polynomial 
basis function. We do not implement this technique here as the goal of this paper is not 
to achieve the most accurate value of bound state energies. However, it can be done in 
principle and could improve the convergence of the numerical method. 

TABLE I. Results of the numerical computation for the ground state of for different mesh sizes 
and B-spline of order 7. Here, are the number of elements in each coordinates while iV* is the 
total number of basis functions utilized. The maximum coordinate was fixed to £ max = 30 a.u. and 
the angular momentum to j z = 1/2. The calculations are to be compared with the results from 



29j where the authors obtained E H + = -1.10264158103 a.u.. 



N 5 


N v 


N* 


E H+ 


(a.u.) 








Min-max 


RR 


8 


8 


182 


-1.102590816884 


-1.102590816895 


10 


10 


240 


-1.102638533873 


-1.102638533934 


12 


12 


306 


-1.102641366239 


-1.102641366228 


14 


14 


380 


-1.102641554428 


-1.102641554501 


16 


16 


462 


-1.102641577089 


-1.102641577085 


18 


18 


552 


-1.102641580210 


-1.102641580229 


20 


20 


650 


-1.102641580782 


-1.102641580825 



B. Spectra of diatomic molecules 

In this section, the spectra of dihydrogen and dithorium are presented. They are shown in 
Tables [TTT1 and HVl for j z = 1/2. The spectra are calculated using a mesh of 14x14 elements 
and 17x17 elements for the min-max method for dihydrogen and dithorium, respectively, 
while a mesh of 30x30 elements is utilized in the Rayleigh-Ritz method. The other param- 
eters are set to the same values as in the last section where the convergence of the ground 
state was discussed. The binding energy values in the mass gap (— mc 2 , mc 2 ), corresponding 
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TABLE II. Results of the numerical computation for the ground state of Th 2 + for different mesh 
sizes and B-spline of order 7. Here, N^ tTI are the number of elements in each coordinates while ./V* is 
the total number of basis functions utilized. The maximum coordinate was fixed to £ max = 15 a.u. 
and the angular momentum to j z = 1/2. The calculations are to be compared with the results from 



29j and 



2J where the authors obtained E Th +_ = -9504.756746922 a.u. and E Th + = -9504.752 



a.u.. 



Nc 




N* 


Erriul79H 

Th 2 


(a.u.) 








Min-max 


RR 


8 


8 


182 


-9503.998584802 


-9504.592867005 


10 


10 


240 


-9504.333585765 


-9504.687658554 


12 


12 


306 


-9504.466070634 


-9504.711111628 


14 


14 


380 


-9504.539502492 


-9504.722791962 


16 


16 


462 


-9504.586247153 


-9504.730034585 


18 


18 


552 


-9504.618392312 


-9504.735005730 


20 


20 


650 


-9504.641636959 


-9504.738611929 


24 


24 


870 


-9504.672557123 


-9504.743429586 


30 


30 


1260 


-9504.698874401 


-9504.747552293 



to bound states, are shifted by mc 2 to have a comparison with non-relativistic results. The 
values in the continua however are not shifted and calculated with the Rayleigh-Ritz method 



only. The results of the dithorium spectrum can be compared to the ones in [37J]. Both are 
generally in good agreement, although a small discrepancy can be seen for the higher excited 
states. 

In the Rayleigh-Ritz method, the n hinding bound state energies shown in Tables UTTl and HVl 
correspond to the 2N+1 to 2A r +l+nbinding eigenvalues of the matrix C (once the eigenvalues 
are ordered in increasing order). The other eigenvalues can be associated to the "discretized" 
negative (the first to the 27V'th eigenvalues) and positive (the 2N + 2 + n binding 'th to the 
4iV'th eigenvalues) energy continua. For the min-max method, the bound state energies 
shown corresponds to the solution of Ak(E) = for the rending lowest energy eigenvalues. 
For the diatomic molecules considered, the spectra calculated with both methods are in very 
good agreement. The small discrepancy remaining is mostly due to the use of different mesh 
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sizes. 

The convergence of the excited states is very similar to one of the ground state: all values 
are approached from above and the order of convergence is close to the one of the ground 
state. The same is true for the states in the positive energy continuum, that is for E > mc 2 . 
For the negative energy states, the convergence occurs from below, but otherwise, follows the 
same trends as the other cases. The energy values in the continua (especially their smallest 
and largest eigenvalues) depend on the size of the domain. In the dithorium calculation, 
the domain was smaller which yielded less accurate value in the continua (not shown in the 
table) but better accuracy of the bound states. In all cases, the eigenvalues of the positive 
and negative energy continua accumulate at the points mc 2 and — mc 2 , respectively. 



C. Spurious states 

The results for the spectra of diatomic molecules presented in the last section showed a 
spurious state in the dithorium spectrum calculated with the naive Rayleigh-Ritz method 
while the other methods did not. Spurious states usually appear as eigenstates with an 
energy in the interval —mc 2 < E < E groun d because of their highly oscillatory behavior. 
This was not observed in the numerical results. Moreover, it was proven mathematically 
that the min-max method is free from these numerical artifacts [18]. The spectra predicted 
by the min-max and the Rayleigh-Ritz methods coincides (up to numerical errors), implying 
that our version of the Rayleigh-Ritz method using kinematically balanced function is also 
free from these unphysical states. 

These last arguments are mainly qualitative. A more convincing approach proceeds by 
computing the spectrum for an atom (by setting Z2 = 0) and by comparing to the well-known 



analytical formula for the atomic binding energy given by 48 



mc 2 , 
E nj = , ~ mc , (53) 

_|_ Z z a' ! 



where n is the principal quantum number, j is the angular momentum and 

5 3 =3 + \-sj(j+ l ^j-Z*a 2 . (54) 

Of course, the numerical methods are not optimized for atomic calculations, but these results, 
albeit not very accurate, allow showing that no spurious states appear. The results for 
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TABLE III. Results of the numerical computations for the spectrum of for a mesh size of 30 x 
30 and B-spline of order 7. The states of the positive and negative continua are computed with 
the Rayleigh-Ritz method and only the first 25 states are shown. 



RminH 

JJUUllU 




Binding energy (a.u.) 




TsJprrati vp 


Positive 




Min-max 


RR 




c f\ n t" l n n n m fan 

LjVJLLUILLU. U.1L1 \ CI . IX . 


continuum (a.u.) 


1 

X 


-1 


1026413662 


-1 


1026415808 


i 

X 


-1 8778 95940 


18778 


86549 


9 


-0 


6675525594 


-0 


6675527718 


9 




18778 


86561 


KJ 


-0 


w a a i—t/t a r r /in 

4287795568 


-0 


4287811584 




-18778 96471 

XOl 1 O. JUt 1 J- 


18778 


86562 


A 


-0 


3608697621 


-0 


O /"* A A l"7 -1 n/inr 

3608710695 


A 


-1 8778 97984 

XOl 1 O.i/i ^Ol 


18778 


86741 


r. 


-0 


2554175614 


-0 


2554197033 


o 


1 ^77^ Q^9^^ 


18778 


86746 


a 
u 


-0 


A « 1 ^7 A A ^7 /"* A A 

2357807609 


-0 


f "V fa ^ 7* 77 A "I A /"* A -| 

2357812681 


6 

U 


-1 S778 QRd7^ 

XOl I O. i/O 1 ! 1 U 


18778 


86808 


7 


-0 


a a r* f—r a a -t A c~\c~\ 

2267021482 


-0 


2267030696 


7 


-1 8778 9Q077 

XOl 1 O.C/JU 1 I 


18778 


86917 


8 
o 


-0 


a a a a /"* a -i ct r* r* 

2008621355 


-0 


2008689095 


o 


-1 8778 9Q390 


18778 


88272 


Q 


-0 


1776816232 


-0 


1776839788 


Q 

U 


1 8778 QQR84 


18778 


88275 


1 n 

1U 


-0 


1373089205 


-0 


-I A ^7 A -| ^ ^7 /~i A /~i 

1373147686 


i n 


1 8778 QQQ7Q 


18778 


88617 


1 1 

X J. 


-0 


1307908409 


-0 


1307928214 


1 1 

X X 


-1 8779 003Q7 

XOl i J.uUOf l 


18778 


88627 


1 9 


-0 


1267066133 


-0 


-| C\r'f—7'\ AA0 1 A 

1267100818 


1 9 

X 


1 877Q 00^44 


18778 


88659 


1 3 
xo 


-0 


-1 A S~* /~l A A ^ -1 

1266438351 


-0 


-1 A /~* J~* A A ~t A A A 

1266441499 


1 3 


-1 8779 00643 

XO 1 I tJ .UUUt:0 


18778 


89067 


14 


-0 


1261986510 


-0 


-i rvA -i a a a A A a 

1261992440 


1 4 


-18779 01959 

XOl li?.UXZ;(Ji? 


18778 


AAA AA 

89068 


1 5 

-LO 


-0 


1158897902 


-0 


1159009024 




1 &77Q 01 ^d9 

-XOl 1 y.UltJ'iZi 


18778 


A A A 

89977 


1 6 
XU 


-0 


1053558675 


-0 


10536H251 


1 fi 


1 &77Q 01 QfH 
-xoi I y.uiyuo 


18778 


A A A A ^7 

89987 


1 7 


-0 


0852450505 


-0 


a a r- - a r - w a a a a 

0852548082 


1 7 


1 £77Q 099^4 


18778 


90406 


1 8 


-0 


0823477309 


-0 


0823523149 


1 8 

-LO 


-1S77Q D9fi^Q 


18778 


90575 


1 Q 

X Cf 


-0 


0804553251 


-0 


0804564514 


1Q 


1 ^77Q 0^400 
-ioi i y.uo^uu 


18778 


90766 


20 


-0 


0802631100 


-0 


0802631662 


20 


-18779.03472 


18778 


91436 


21 


-0 


0802102415 


-0 


0802110614 


21 


-18779.03980 


18778 


91438 


22 


-0 


0802047983 


-0 


0802048234 


22 


-18779.04047 


18778 


91544 


23 


-0 


0800201297 


-0 


0800252078 


23 


-18779.04625 


18778 


91591 


24 


-0 


0730502242 


-0 


0730676123 


24 


-18779.04822 


18778 


91599 


25 


-0 


0649840057 


-0 


0649993038 


25 


-18779.05030 


18778 


91605 
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TABLE IV. Results of the numerical computation for the spectrum of Th 2 . The mesh size is 
indicated on the second line. The B-splines are of order 7. 



O LcXlA. -ft 


Naivp RR 




RR 


l\/T m_ m a "v 

1VX1J.1X1J.CLA. 




14 x 14 


30 x 30 


30 x 30 


16 x 16 


1 

_L 
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Q^ftd 794^99^ 


Q^ifM 747 c W : i 
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41 97 87QQ c i' : U 
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R 

u 
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94^ Q^^Q 1 ^ 
-i^OO.JjO / voo 
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7 
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-901 4391 1 03 
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o 


1 Q1 R ^97^474 
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1 Q1 ^ fi7fi1 9fi7 


q 
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1 fUQ 9Q9Q1 48 




1 fi43 Q390fifi^ 


1 


- X O L ± O . O O £ a U O f± 


1 344 08^870 
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1 31 3 8071 Q1 (\ 


1 31 3 7fif)fi8QQ 


ii 


-1339.1123032 


-1333.5368147 


-1303.6850950 


-1303.6580541 


spurious 


-1218.2113620 


-1204.6990945 






12 


-1169.3956263 


-1159.1761393 


-1089.6415827 


-1089.6356220 


13 


-1138.5709512 


-1131.0151665 


-1084.3699127 


-1084.3519981 


11 


-1046.2053120 


-1045.4764538 


-1028.1920826 


-1028.1912423 


15 


-1018.4013912 


-984.5252901 


-969.6816867 


-969.64172165 



the spectrum of Th 89+ are shown in Table |Vj For all eigenenergies considered, there is 
always a one-to-one correspondence between the analytical and the Rayleigh-Ritz results, in 
contradistinction with the spectrum obtained from the naive Rayleigh-Ritz method which 
exhibits spurious states. 
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TABLE V: Results for the spectrum of Th 89+ . The mesh size is 30x30 and the B-splines are of 
order 7. The states are denoted in spectroscopic notation. 
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2 
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5gs 
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IV. CONCLUSION 

In this work, we presented two numerical methods to solve the single particle time- 
independent Dirac equation. The first one was based on a min-max variational principle 
while the second one used a combination of the Rayleigh-Ritz method and kinematically 
balanced basis functions. For comparison purposes, we also included a description of the 
naive Rayleigh-Ritz method. All were based on a B-spline basis function discretization which 
allowed obtaining a high accuracy and a sparse matrix structure in the discretized equations. 
We applied these methods to the computation of the two-centre Coulomb problem ground 
state energy and spectrum. Because of its axial symmetry and simple structure, it was 
convenient to use prolate spheroidal coordinates. These techniques were used specifically 
to compute the spectra of the molecule and the quasi-molecule Tli2 79+ . A comparison 
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with results in the literature for the ground state demonstrated that our methods yield very 
accurate and convergent results, especially for dihydrogen. More importantly, no spurious 
states were reported in these numerical schemes and thus, they both could be used to evaluate 
radiative QED corrections for diatomic molecules which necessitate sums over intermediate 
states. This conclusion was reached by comparing the calculated spectra of dithorium and 
thorium to results obtained from the naive Rayleigh-Ritz method. 

The two methods have strengths and weaknesses. In terms of computation time, the 
Rayleigh-Ritz methods were much faster, especially for the computation of the whole spec- 
trum. This happens because in the min-max method, the solution of the eigenvalue equation 
Afc(-E') = necessitates many iterations (typically between 20 and 30) to obtain a decent 
accuracy and each iteration requires a solution of the eigenvalue problem. This could be 
improved somewhat by using an iterative eigensolver optimized for the computation of few 
eigenvalues. Therefore, if one is only interested in the computation of the first few excited 
states while using a very large mesh, it may be advantageous to use the min-max method 
combined with a version of these iterative eigensolvers. In terms of accuracy, both methods 
yielded very similar results, although the convergence was slightly better for the Rayleigh- 
Ritz method in the dithorium case. This is in contradiction with the conclusion reached 
where the min-max method showed much better accuracy. This discrepancy may 



m 



be explained by a slightly different choice of basis functions (see ( 133|) versus (6) of [381]). 
Nevertheless, the main advantage of the Rayleigh-Ritz method is the fact that it gives the 
whole spectrum directly from the solution of the generalized eigenvalue problem. 

These methods can then be utilized in many applications. Among others, this can be 
used to investigate relativistic laser-matter interaction: the solution obtained from these 
methods can be used as an initial condition for the solution of the time-dependent Dirac 
equation. This will be the subject of future investigations. 



Appendix A: Explicit expression for min-max method 



The explicit expression of matrices An, Aw and A 12 is obtained by using the Dirac 
equation in cylindrical coordinates given in 41] . Then, combining the ansatz in Eq. f fT0|) 
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with the basis function expansion, we get 



[An].- = / 



d 3 x 



(d z B^)(d z Bf) + (d r Bn(d r B^) 



[A 22 ] 



[a 12 ; 



^l/Q n(l)^(l) 



£ + mc 2 — V 



d 3 x 



+(V + mc 2 — E)Bf'B { j 



(d z BP)(d z Bf) + (d r B^)(d r Bf) 



+ ^B^Bf + ^\d r Bf) 
+ ^{d r B? ) ) B f 



d 3 x 



E + mc 2 -V 
+ (V + mc 2 - E)Bf ) Bf 



(d z B^)(d r Bf) + ^\d z Bf) 



-(d r B^)(d z Bf) + ^(d z B^)Bf 

^2 



X ■ 



E + mc 2 -V 



(Al) 



(A2) 



(A3) 



The last expression can then be expressed in prolate spheroidal coordinates by using 



and the integration measure is given by 



(A4) 
(A5) 



d 3 x = R 3 (£ 2 - r] 2 )d^dr]d9. 



(A6) 



Appendix B: Explicit expression for the naive Rayleigh-Ritz method 



The explicit expression of matrices C^, and is obtained by starting with the 
Dirac equation in cylindrical coordinates. By assuming that the basis functions are the same 
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for the large and small components, we get 



We also have 



[c& 

[^22 

[C& 

[^22 

[Cg 

[C 

[C (3 

[c 
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J d 3 x^(V + mc 2 )B^Bf^, 



d\r{ {V + mc 2 )Bl 2) Bf 
(fxliV-mc^B^B^ 



d 3 x{ {V -mc 2 )Bf ) B J 



d 3 xlBi l) d z Bf 



- I d 3 xi Bf ] d z Bf 



j d 3 x< 




d r - 


/ii 

r 




1 d 3 x< 




d r + 


r 


Bf] 



(Bl) 

(B2) 

(B3) 

(B4) 

(B5) 

(B6) 

(B7) 

(B8) 
(B9) 



PS] 



d 3 x{ B - B 



(i) 



d 3 x{ BY B 



3 



(2) o(2) 



(BIO) 

(Bll) 
(B12) 



The last expressions can then be expressed in prolate spheroidal coordinates with (lA4j) .( lA5]) 
and (jMp . 



Appendix C: Explicit expression for the Rayleigh-Ritz method with KBBF 



The explicit expression of matrices An, A 22 and A i2 is again obtained by using the Dirac 



equation in cylindrical coordinates given in 



41] . Then, combining the ansatz in Eq. f fTUj) 
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with the basis function expansion, we get 



[ C 22] ij 



[C22] ij 



d 3 x^(V + mc 2 )B ( i 1) Bf 
d s x^(V + mc 2 )B^Bf 

rf 3 x| (d M B^)(d z Bf) + (drB^drB^) 



-^(d r B^)Bf 



(V — mc 2 ) 
4m 2 c 2 



(fx 



(d x B^)(d z Bf) + (drB^)(drBf) 



+ ^ B (2) B (2) + ^ B (2) {drB (2 }] 



r 



+ ^(d r B^)Bf 



(V — mc 2 ) 



d 3 x 



4m 2 c 2 

(dM 1] Kd r Bf) + ^\d z Bf) 



r 

-(d r B^)(d z Bf) + !±{d z B?>)Bf> 



x 



d 3 x 



(y - mc 2 ) 
4m 2 c 2 

(8,b, <1) )(8,b}") + (arrows]": 



+^Bf ) Bf + ^Bf\d r Bf ) ) 

[C? 2 % = J d*x{(d z BU)(d r B?) + ) 
-(^)M 2 >) + ^(^)S«' 
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We also have [C®].. = [Cg] .., [Sg% = 2m[Cg J ] .., [Sg% = 2m[Cg] y and [Sg] 



,(3)1 



s(3) 



(3)1 



J (3)1 



,(3)1 



:(3)1 



2m [Cg]... The last expressions can then be expressed in prolate spheroidal coordinates 
with (|A4jh(|A5|) and 
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